Geographic data in Python

What is geographic data?

• Thus far: we’ve been working with data (largely) in tabular data frames, where each column represents data “attributes”
• Geographic data: includes information on location - specifically, where the observation is located on the Earth’s surface

Maps!

Longitude and latitude: the basics

Geographic coordinates

• Longitude (X) and latitude (Y) coordinates can be expressed as:
  • Degrees minutes seconds (e.g. \(97^{\circ}21'37''W\), \(32^{\circ}42'38''N\))
  • Decimal degrees (e.g. \(-97.36\), \(32.71\))

Conversion between DMS and decimal degrees:

\[DD = D + \frac{M}{60} + \frac{S}{3600}\]

Coordinate systems

• Geographic coordinate system: coordinate system based on latitude and longitude (coordinates on a sphere)
• Projected coordinate system: coordinate system projected onto a two-dimensional surface (coordinates on a plane)
• Example: http://bl.ocks.org/mbostock/raw/5731632

Map projections

Tiled mapping and Web Mercator

Map projections and distortion

Examples:

• http://bl.ocks.org/curran/raw/01aa2685f083b6c1b9fb
• http://bl.ocks.org/enjalot/raw/bd552e711b8325c64729

Vector data

Raster data

Vector data: CSV files

• Example: Starbucks locations in Chicago

import pandas as pd

starbucks = pd.read_csv("http://personal.tcu.edu/kylewalker/data/chicago_starbucks.csv")

starbucks.head()

GeoPandas

• GeoPandas: Python package that extends pandas to represent geographic objects as GeoDataFrames

• Key column: geometry, which represents the geometry type of the data (point, line, or polygon) and the sequence of coordinates

Making a GeoDataFrame

import geopandas as gp

sbgeo = gp.GeoDataFrame(starbucks,
                        geometry = gp.points_from_xy(starbucks.Longitude,
                                                     starbucks.Latitude),
                        crs = 4326)

Making a GeoDataFrame

• gp.GeoDataFrame(): function used to generate the GeoDataFrame, which can take an existing pandas DataFrame as its first argument

• Parameters:
  • geometry: how to represent the data as a geographic object. We use the gp.points_from_xy() function here to “make” the geometry from existing Longitude and Latitude columns.
  • crs: a code that specifies the dataset’s coordinate reference system. Most coordinate systems can be represented by 4- or 5-digit codes (see https://spatialreference.org/)

Mapping point data

sbgeo.plot()

Adding a “basemap”

import contextily as ctx

# "Project" to Web Mercator, used by tiled mapping services
sbgeo2 = sbgeo.to_crs(3857)

# Use the basemap function found at https://share.cocalc.com/share/56cef78b3fda5e0338c2da71b222e6b275a9193e/contextily.ipynb?viewer=share
def add_basemap(ax, zoom, url='http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'):
    xmin, xmax, ymin, ymax = ax.axis()
    basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=zoom, url=url)
    ax.imshow(basemap, extent=extent, interpolation='bilinear')
    # restore original x/y limits
    ax.axis((xmin, xmax, ymin, ymax))

Adding a “basemap”

p = sbgeo2.plot()
cx.add_basemap(p, zoom = 11)

Vector data: GeoJSON

• Extension to the JSON format that encodes geographic coordinates for datasets

• Mexico’s states as GeoJSON

Mapping polygons

import geopandas as gp

mxgeo = gp.read_file("https://gist.githubusercontent.com/walkerke/76cb8cc5f949432f9555/raw/363c297ce82a4dcb9bdf003d82aa4f64bc695cf1/mx.geojson")


sns.set(style = "white")
mxgeo.plot(column = "pri10", cmap = "RdPu", figsize = (10, 8))

Mapping polygons

Should you map?

• Beware the “amazing map”…

Maps vs. charts

• For discussion: which visualization is more effective for showing differences in our data?

Chart for comparison

Code to reproduce the chart

import pandas as pd, seaborn as sns, matplotlib.pyplot as plt
sns.set(style="whitegrid", font_scale = 1.3)
df = pd.read_csv("http://personal.tcu.edu/kylewalker/mexico.csv")

plt.figure(figsize = (10, 8))

p = sns.stripplot(data = df.sort_values('pri10', ascending = False), 
                  x = 'pri10', y = 'name', palette = "RdPu_d", 
                  orient = 'h', size = 8)
p.set(xlabel = "% of workforce in primary sector", 
      xlim = (0, 50), ylabel = "")
p.axes.xaxis.grid(False)
p.axes.yaxis.grid(True)
sns.despine(left = True, bottom = True)

Raster data in Python

• To work with raster data in Python, use the rasterio package

import rasterio as rio

# Dataset of elevation in the Tarrant County area
elev = rio.open("http://personal.tcu.edu/kylewalker/data/tarrant_elevation.tif")

Plotting raster data

from rasterio.plot import show, show_hist
sns.set(style = "white")
show(elev)

Raster data: histograms

show_hist(elev, bins = 100)